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We present a probabilistic theory of random walks in turbid media with non-scattering regions. 
It is shown that important characteristics such as diffusion constants, average step lengths, cross- 
ing statistics and void spacings can be analytically predicted. The theory is validated using Monte 
Carlo simulations of light transport in heterogeneous systems in the form of random sphere packings, 
and good agreement is found. The role of step correlations is discussed, and differences between 
unbounded and bounded systems are investigated. Our results are relevant to the optics of hetero- 
geneous systems in general, and represent an important step forward in the understanding of media 
with strong (fractal) heterogeneity in particular. 



I. INTRODUCTION 

Multiple scattering of light is ubiquitous in nature, and 
a proper account of the phenomenon is essential in impor- 
tant areas such as atmospheric science [IH2], biomedical 
optics [4]-[6] and material characterization [TUTU] . Despite 
the abundance of situations in which scatterers are not 
uniformly distributed, multiple scattering of light is of- 
ten treated within the framework of radiative transfer 
and diffusion theory while assuming exponentially dis- 
tributed ballistic segments. At the same time, significant 
efforts have also been directed towards the understand- 
ing of heterogeneous systems and non-Poissonian trans- 
port mechanisms pTUT4] . In fact, the presence of hetero- 
geneities has important implications in a wide variety of 
contexts, from the optics of cloud systems [3j H5HT7] to 
spectroscopy of the human body [18-26 , to transport in 
colloids and foams [27-30 . However, since light trans- 
port in such cases is intrinsically related to the specific 
heterogeneity, its description in general terms remains 
elusive. 

Here, we present a probabilistic theory for light trans- 
port in turbid media with non-scattering regions in- 
side. We show that important characteristics such as 
asymptotic diffusion constants, step distributions and 
void crossing statistics can be analytically predicted via 
simple formulas. The success of the theoretical approach 
is verified by analysis of random walks in random three- 
dimensional polydispersive sphere packings (the inter- 
sphere volume set to be homogeneously turbid). Our 
findings provide important insight on light transport in 
materials with spatial heterogeneity in scatterer density, 
including complex porous materials and disordered opti- 
cal materials with engineered fractal heterogeneity such 
as Levy glasses 



II. ANALYTICAL DESCRIPTION OF HOLEY 
RANDOM WALKS 

Let us start to consider the general case of a ran- 
dom walk in a composite media consisting of arbitrar- 
ily shaped non-scattering regions embedded in a turbid 
medium (a holey system), the void filling fraction being 
<f> (cf. Fig. [I). 
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FIG. 1. A holey random walk. Random walkers change direc- 
tion randomly in time, but only when traveling in the turbid 
medium found in between the non-scattering regions (holes). 
Transport is governed by the resulting step length distribution 
and step correlations in a complex manner. 

We will assume that isotropic scatterers are randomly 
distributed in the turbid medium, meaning that the dis- 
tance between isotropic scattering events is exponentially 
distributed with a scattering mean free path l s . The 
scattering mean free path is related to single scattering 
properties via 
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where n denotes the number density of scatterers and 
<j s the scattering cross section. In general, the random 
step length S between two subsequent scattering events 
in the composite heterogeneous medium consists of two 
parts: a length 5t U rbid £ Exp(£ s ) inside the turbid inter- 
void medium, and a length ^Void inside non-scattering 
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void(s): 



S = Si 



turbid 



s. 



void- 



(2) 



Here, it should be noted that Sturbid and SVoid are n °t 
independent random variables. On the contrary, they are 
positively correlated: a long step in the turbid medium is 
more likely to be accompanied by one or several crossings 
of non-scattering regions. 



A. Equilibrium considerations 

Assuming that the density of states is homogenous 
(which in optics corresponds to a constant refractive in- 
dex), random walkers at equilibrium should sample the 
two constituents according to their respective volume 
fraction. We then expect that a fraction </> of the average 
step is located to voids, i.e. 



E[S mid ] = <t>E[S] 
E[S turbid ] = (1 - <f>)E[S] 

In terms of expectations values, we have 

E[S] = E[S vo id] + ^[Sturbid], 



(3) 
(4) 



(5) 



and using that we per definition also know that 
^[Sturbid] = £ s , we reach E[S] = (j)E[S) + £ s and con- 
clude that 



E[S] 
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This important result that tells us two things: the mean 
step length is independent of the shape and size distri- 
bution of the non-scattering regions, and it equals the 
homogenized scattering mean free path £^ that would 
govern a system where the scatterers were randomly dis- 
tributed over the full volume instead of only a volume 
fraction 1 — (n going from no to no (1 — </>) in Eq. [I]). 
E[S] is thus dependent only on 0, not on the sizes of the 
voids. This means that the mean step will only grow with 
an increase of 0, and not with the inclusion of larger and 
larger voids. 



B. A general remark regarding diffusion 

Although the mean step is a very important quantity, 
transport properties are determined by the step length 
distribution as a whole. In particular, when steps can 
be considered independent, the diffusion constant that 
governs the macroscopic spreading of random walkers is 
determined by the ratio of the first two moments of the 
step distribution (when finite). This important relation 
can easily be derived by considering the summation of in- 
dependent d-dimensional isotropic random steps, all fol- 
lowing the same (arbitrary) step length distribution. One 



finds that the mean square displacement asymptotically 
grows as 2dDt with a diffusion constant D given by 



D 



E[S 2 ] 
2d XVX E[S] 



(7) 



where S denotes the scalar length (norm) of the individ- 
ual steps, and v is the walk velocity. For purely expo- 
nential steps in 3D, corresponding to a homogenous dis- 
tribution of scatterers, this expression reduces into the 
famous relation 



D = -v£ t 
3 



(8) 



where £ t is the so called transport mean free path (the 
mean length of the exponential steps). It should be noted 
that this relation is normally derived via the radiative 
transfer equation [37j [38] , not by considering summation 
of random variables (but see, e.g., [39l [40] for similar 
considerations). For a general holey system, with het- 
erogeneous distribution of scatterers, steps are not expo- 
nential distributed and the simple formula of Eq. [8] does 
not hold. Interestingly, we have already shown that the 
mean step length, E[S], of the holey system is indepen- 
dent of the scatterer distribution. The second moment, 
E[S 2 ], will on the other hand depend on the particular 
heterogeneity and scatterer density in a complex manner. 



C. Sphere packings as holey systems 

As indicated above, transport in a general holey system 
is complicated. Beside the issue of step length distribu- 
tion, step correlations may be important and must not 
be forgotten. To gain insight on this complicated mat- 
ter we will turn to systems where the non-scattering part 
of the system have the form of a polydispersive random 
sphere packing. In this and the following section, we will 
present an analytical theory of transport in such systems. 
In subsequent sections, this theory will be compared to 
Monte Carlo simulations of random walks in sphere pack- 
ing realizations. 

Let Ti (i = 1, . . . , M) denote the different radii of the 
involved spheres, and n^ their respective number density. 
A random walker traveling in this system will scatter ran- 
domly in time, but only when being outside of the non- 
scattering regions. When being scattered, the chance of 
crossing at least one sphere in the coming step depends 
on the probability that the next step is longer than the 
distance to the next sphere surface (along the new walker 
direction). This distance can be seen as a random vari- 
able, and is here denoted A ps (ps as in point-sphere). If 
a sphere is crossed, the probability to cross also a sec- 
ond sphere before being scattered now depends on the 
sphere-sphere spacing along the direction of the walker. 
Also this distance is a random variable, and we denote 
it A ss (ss as in sphere-sphere). The actual step length 
distribution will be determined by the distributions of 
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these random variables (of course, in combination with 
the inter-sphere scattering mean free path £ s ). Fig. [2] 
illustrates the definition of these two essential distance 
variables. 




FIG. 2. From a random walk perspective, two central charac- 
teristics of the sphere packing are the distribution of (i) the 
distance from a random point r p (scattering event) along a 
random direction u to a sphere surface, A ps (r p ,u), and (ii) 
distance to the next sphere along a random direction when 
leaving a sphere at a random point r s , A ss (r s , u) . Their dis- 
tributions can be assessed via statistical analysis of sphere 
packings. Later, we show that the mean value of A ss , as 
experienced by the random walker, can be analytically pre- 
dicted. 

The average spacing between spheres A ss , as experi- 
enced by the random walker, can be analytically calcu- 
lated from the filling fraction and the sphere distribution. 
If we consider a random line drawn through the sphere 
packing, it is evident that a fraction 1 — of it will be 
drawn through the intersphere medium. Letting E[(] de- 
note the average length of the individual chords through 
the spheres, we have that 



E[A S 



E[A SS ] + E[(] 
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and, similarly to Olson et al. |H], conclude that 



E[A S 



-E[(}- 



(9) 



(10) 



The average chord length E[Q can, in turn, be calculated 
from the sphere distribution. Let pi denote the proba- 
bility that a random chord is made in a sphere of radius 
T{. Intuitively, this probability should be proportional to 
the total surface area of the sphere category, i.e. 



Pi = 



j 



2 

x rf 



(11) 



In fact, this view is equivalent to the equilibrium argu- 
ment adhered to in this article: the total path made 
through the different sphere categories should equal their 
fraction of the total volume. Assuming an isotropic flux 



of random walkers close to the sphere surfaces (see dis- 
cussion in [42 ) , the length of a chord through a sphere of 
radius follows a triangular probability density function 
fi(x) = xj{2rf) in < x < 2r^. The average chord being 
ti = 4r^/3 and the mean square chord being 2rf. The 
probability density function of a random chord ( in the 
polydispersive packing, f^(x), is a weighted sum of the 
individual chord distributions ff. 



Pi x fi(x) = ^2 



Pi x 



i:2ri >x 



2r?' 



The mean overall chord then becomes 

and its mean square 

£[C 2 ] = V> x 2r?. 



4r%- 



(12) 



(13) 



(14) 



D. The exponential spacing model 

From statistical analysis of random sphere packings, 
we generally find that i?[A ps ] « £^[A SS ] and that distribu- 
tions have near-exponential decay. We therefore propose 
a model in which the actual (complicated) distributions 
of A ps and A ss are modeled with a single exponentially 
distributed random spacing A with mean spacing A, i.e. 
A G Exp(A). In fact, it has been shown that spacing 
between random objects is, to a very good approxima- 
tion, exponential at low filling fractions [41]. At higher 
filling fractions, spacing distributions are no longer be 
perfectly exponential, but rather a complicated function 
of the exact structure [41, 43, 44 j. But also in such cases, 
we will show that the errors due to the use of a simple 
exponential spacing model can be reasonably small. 

Within the exponential spacing model, the probability 
that an exponentially distributed step 5t U rbid £ Exp(/ S ) 
will not take the random walker to a sphere is 



Po = P(S tUYhid < A) 



A + 4 



(15) 



Given the memoryless property of the Poisson process 
that governs scattering events, the number of crossings 
naturally follows a geometric distribution. This means 
that the probability to cross n spheres in between two 
scattering events is 



p n = a - p ) n p . 



(16) 



The average number of sphere crossings per step N is, 
within this model, 



E[N] = J2^xP n 



1 



A 



(17) 



In order to fulfill the equilibrium condition, A must be 
chosen so that the average step samples the medium ac- 
cording to volume fractions. Assuming that the number 
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of sphere crossings and involved chord lengths are inde- 
pendent, we should have that 



E[N}E[C] = , 
E[N]E[(]+£ S 9 

Using Eq. [TTJ we then reach the conclusion that 

1 



E[A S 



-E[(] 



(18) 



(19) 



The exponential spacing model we have now outlined 
allows us to the estimate step length distribution. In 
particular, we can use it to estimate the mean squared 
step E[S 2 ]. Again, we will assume that, for each step, the 
involved chord lengths are independent of the number of 
sphere crossings. By splitting the outcome of S 2 into 
the number of crossings and summing the conditional 
expectations (details in Appendix [A]) , we find that 

E[S 2 } = E[S 2 urhid ] + E[S 2 oid ] + 2£;[S t urbid£void] 
= 2£ 2 



J2PnX(nE[( 2 ]+n(n-l)E[(] 2 ) 



n=0 



2^P„xn£[(] x (n + 1) 
E[e\ 



i=0 



A + £ s 



= 2ii + 04 x 



E[(] 



(20) 



E. The diffusion constant for holey system 

Assuming that step correlations are negligible, and re- 
membering that E[S] = £h, the expression for E[S 2 ] 
given above (Eq. |20| ) allows us to write a closed form 
expression for the diffusion constant that governs macro- 
scopic transport. Introducing the homogenized diffusion 
constant 



1 



(21) 



and the "chord-domain" diffusion constant (the diffusion 
constant for a random walk with steps purely following 
the chord step distribution) 



D 



c 



1 E[C 2 } 



(22) 



we reach, via Eq. [7j the surprisingly simple relation 

D = D h + 0D C (23) 

Interestingly, the diffusion constant of the heterogeneous 
system is given by the homogenized diffusion constant 
plus a term which is independent of the inter-void scat- 
tering mean free path £ s . Note also that, as expected, the 
expression reduces to the homogenized diffusion constant 
as (j) approaches zero. 



III. MONTE CARLO SIMULATIONS OF 
TRANSPORT IN QUENCHED DISORDER 

To illustrate the validity and importance of the the- 
ory presented in earlier sections, this section will com- 
pare it with Monte Carlo (MC) simulations of transport 
in quenched random sphere packings. We use the term 
quenched to emphasize that the random walk is done in 
the quenched (frozen) heterogeneity of a random sphere 
packing, and not in fully annealed disorder model. The 
sphere packings used in our simulations are created by 
random sequential addition of spheres [45] in descend- 
ing order of size. Moreover, they are created to have 
periodic boundary conditions, i.e. so that the complete 
(cubic) sphere packing can act as a unit cell that can be 
stacked in all directions. As illustrated in Fig. [3j random 
walks are performed so that when crossing a boundary 
of the unit cell, the random walk can be continue on the 
opposite side of the same unit cell (while keeping track 
of unit cell coordinates). The use of periodic boundary 




FIG. 3. 2D illustration of unbounded holey random walks 
based on a sphere packings with periodic boundary condi- 
tions. When a walker leaves the sphere packing unit cell 
(marked square), the walk continues on the opposite side of 
the same sphere packing. The path of interest (red) is re- 
constructed by keeping track of boundary crossings, but path 
outside the unit cell is in reality a path inside the one and 
only unit cell (blue path) . The use of periodic boundary con- 
ditions in this manner enables the use of realistically sized 
sphere packings, efficient averaging over disorder realization 
by allowing walkers to start also close to the sphere packing 
boundaries, and long-time tracking of spreading. At the same 
time, the size of the unit cell is made large enough to make 
effects of periodicity negligible. 

conditions in this manner allows us to use realistic sizes 
of sphere packings, and still being able to track walkers 
for long times and average dynamics over local disorder 
realization in an efficient way. Still, the size of the unit 
cell is made large enough to render effects of periodic- 
ity negligible (sizes of utilized sphere packings are care- 
fully stated). We base our investigations on equilibrium 
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starting conditions, i.e. allowing walkers to start ran- 
domly over the whole unit cell (including the interior of 
the non-scattering regions). For step statistics, only full 
steps between scattering events are considered. To be on 
the optical time scale, walkers are set to travel at the 
speed of v = 200 um/ps (corresponding to a refractive 
index of 1.5). 



A. A monodisperse sphere packing 

The first case study is a holey system based on a 
monodispersive sphere packing. The non-scattering part 
of the material is constituted by randomly placed spheres 
of radius 50 |im, the sphere filling fraction being <p — 0.3. 
Fig. [4] exemplifies the mean square displacement ob- 
tained by averaging the dynamics of 5 x 10 5 random walk- 
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FIG. 4. Dynamics for a holey system based on a monodisper- 
sive sphere packing (r = 50 (xm, <j) — 0.3, £ s = 200 \im and 
v = 200 pm/ps). The diffusion constant extracted from anal- 
ysis of the MSD evolution is in excellent agreement with the 
diffusion constant predicted by our theory (Eq. 23). For this 
system theory predicts a diffusion constant of 19798 |im 2 /ps, 
and linear fit over the temporal range indicated in gray in 
part (a) gave 19780 =b 34 [xm 2 /ps (mean ± standard error of 
the mean, as obtained from five simulation repetitions with 
10 5 random walkers each). The middle graph (b) shows a 
log-log plot of (a). By showing the slope of the log- log plot, 
the bottom graph (c) reveals the transition from ballistic to 
diffusive dynamics and verifies that the temporal range used 
in our analysis falls within the asymptotic diffusive limit. 



ers launched randomly in a cubic system with a side of 
about 12 mm (1025017 randomly placed spheres). The 
average spacing between spheres is, in this case, 155.6 |im 
(as given by Eq. 10). As discussed in the figure caption, 
the outcome is in excellent agreement with our theory. 

The applicability of our theory is further elaborated 
in Fig. [5j in which we study how the diffusion con- 
stant depend on the inter-sphere scattering mean free 
path £ s . Interestingly, the example shows that even 
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FIG. 5. The analytical theory for the diffusion constant (red 
solid line, Eq. 23) is in excellent agreement with quenched 



Monte Carlo simulation of transport in a monodispersive 
sphere packing with r = 50 U-m and <j> — 0.3 (solid dots, 
average of 5 repetitions of 10 5 walkers per level of £ s ). This 
graph shows how D obtained from simulations exceed the 
homogenized diffusion constant Dh (dashed line) exactly by 
the amount <j)D^ , largely independent of £ s . As scatter- 
ing becomes stronger and stronger (reduction in £ s ), anti- 
correlations between steps should reduce D, and deviations 
from our theory is expected. The bottom graph clarifies that, 
in this particular case, this effect becomes important only for 
£ s smaller than approximately r/2 = 25 U-m (here, all five 
repetitions are shown). Note that the mean spacing between 
spheres is .E[A SS ] « 156 fim. 

for cases where the average spacing between spheres, 
2£[A SS ], is several times larger than £ S: step correla- 
tions are negligible and our analytical theory is appli- 
cable. Of course, when i s becomes significantly smaller 
than the void size, the increase in return probabilities 
should affect transport (steps become anti-correlated) . 
This in good agreement with the result that step corre- 
lations in this monodispersive system is important only 
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for £ s < r/2 = 25 |im. The strong quenching regime 
requires separate treatment and is beyond the scope of 
this article. 



B. A fractal sphere packing 

Let us now turn to the more general case of a polydis- 
persive sphere packing, i.e. a situation where the non- 
scattering regions vary in size. We will study the limit- 
ing case of a system with fractal heterogeneity (over two 
orders of magnitude) and a high sphere filling fraction. 
More specifically, our sphere packing is based on M = 18 
sphere categories where the radii T{ are sampled exponen- 
tially sampled from r m i n = 2.5 |Ltm up to r max = 200 (am, 
i.e. 



Probability density 
0.2 



x exp log 



1 



M - 1 



(24) 



The number density of the spheres is set proportional 
to 1/rf (i.e., the different sphere categories all occupy 
the same volume fraction), and the filling fraction of the 
system as a whole is <\> = 0.7005 (9507676 spheres in a 
cube with a side of about 2 mm). 

Fig. [6] highlights the success of the exponential spacing 
model also for such a complex heterogeneous system. In 
particular, this means that we have good quantitative 
knowledge on the step length distribution (including the 
sphere crossing statistics). 

Using i s = A as a first test case, Fig. [7] shows the 
onset of diffusion and the asymptotic diffusion constant. 
Interestingly, the asymptotic diffusion constant is in this 
case about 25% lower than what we would have if steps 
were independent. With decreasing scattering strength, 
the importance of step correlations steadily decrease, and 
the diffusion constant becomes in good agreement with 
our analytical expression D = Dh + <j)D^ . This behavior 
is illustrated in Fig. |8j which shows simulations for £ s be- 
ing 1 to 8 times A. In this context, we want to emphasize 
that step correlations per se do not affect the applicabil- 
ity of the exponential spacing model and the accuracy 
of the related E[S 2 ) estimation (Eq. 20). Instead, they 
only affect the validity of the assumption of independent 
steps that leads to D = D h + </>L> c (Eq. [23] ). 

Before ending this section, we would like to make a 
brief comparison to the monodispersive system studied 
in the previous section. There, we found that step anti- 
correlations started to play a role only when the scat- 
tering mean free path was smaller than the size of the 
involved spheres. For our polydispersive system, the situ- 
ation is more complicated. The average spacing between 
spheres are only a few |im, and the sphere sizes range 
from 5 to 400 |im. One can expect that large spheres 
will give rise to strong step anti-correlations, but since the 
transport between these larger spheres is a complicated 
function also of the smaller spheres in-between, it is diffi- 
cult to know for which t s this effect becomes important. 
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FIG. 6. The theoretical predictions of the exponential spac- 
ing model is in good agreement with both quenched simu- 
lations and statistical analysis of a strongly polydispersive 
sphere packing. The success is here illustrated for the case of 
a fractal sphere packing with M — 18 radius categories rang- 
ing from r m in = 2.5 |^m to r max = 200 \±m (9507676 spheres 
in a cube with a side of about 2 mm, = 0.7005). As shown 
in the top graph, the analytically predicted mean spacing of 
A = 5.89 |^m is in good agreement with spacing distribu- 
tions obtained via statistical analysis of the sphere packing 
CE[A ps ] « E[A SS ] « 6 pm). The bottom graph highlights 
that the crossing probabilities P n predicted by the exponen- 
tial spacing model (Eq. 16 ) agree almost perfectly with sim- 
ulations running at £ s = A. In particular, setting t s = A 



(as done in this simulation) clearly renders E[N] = 1. The 
mean and mean squared steps are also in good agreement 
with theory. Theory predicts E[S] — £h — 19.669 |^m and 
E[S 2 ] = 1776.7 pm 2 , and simulations gave 19.679 ±0.004 \im 
and 1692±1 pm 2 , respectively (mean ± standard error). That 
the observed mean squared step is a few percent lower than 
the theoretical prediction is consistent with the fact that the 
exponential model slightly overestimates the number of higher 
order multi-crossings (n > 3). This, in turn, can be un- 
derstood by looking in the upper graph and noting that the 
sphere-sphere spacing distance often will be underestimated 
by the exponential spacing model. 



As shown in Fig. [8j step correlation is of little important 
once i s is larger than, say, 4A « 25 |Ltm. As mentioned 
also in the figure caption, this scattering strength gives 
a mean step E[S] = £h of about 80 |Ltm. That correla- 
tions are weak can then be (partly) understood by realiz- 
ing that most spheres are significantly smaller than this 
mean step length (the mean chord of the system is, in 
fact, not longer than E[(] = 13.8 |Ltm). 
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FIG. 7. Transport dynamics for a holey system based on a 
fractal sphere packing (sphere radii ranging from 2.5 |^m to 
200 \im, each category occupying about 4% of the volume, the 
overall sphere filling fraction being (f) — 0.7005). Here, the 
walker crosses on average one sphere in-between scattering 
events (£ s — A = 5.89 urn). The exact E[S 2 ] as obtained 
from simulations would, if steps were independent, result in a 
diffusion constant of 2866 uxn 2 /ps (step data presented in Fig. 
[6] but note that exponential spacing model in this particular 
case overestimates the E[S 2 ]/E[S] ratio by about 5%). In 
contrast, a linear fit over the temporal range indicated in gray 
in part (a) gives D = 21 11 ±7 u.m 2 /ps (mean ± standard error 
of the mean, as obtained from five simulation repetitions with 
10 5 random walkers each). That the observed D is about 
25% lower is clear evidence that step correlations play an 
important role. 
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FIG. 8. Also for strongly poly dispersive systems, our prob- 
abilitisic theory is successful. The role of step correlations 
is, however, a complicated matter. In our particular frac- 
tal system, where 400 \im diameter spheres are the largest 
non-scattering regions, the role of correlation becomes small 
when i s is a few times larger than the mean void spacing 
A. At i s = 4 A « 24 urn, the asymptotic D obtained from 
MSD evolution is then only 4% lower than what the step 
statistics from quenched simulations tells us. The mean step 
E[S] — in is then about 80 \im, a length which apparently is 
large enough to render step anti-correlation related to cross- 
ing of large spheres relatively small. Here, it should be noted 
that the exponential spacing model overestimates E[S 2 ] by a 
few percent (cf. discussion in Fig. |6|. The exact diffusion 
constant we would see if steps were independent can be esti- 
mated from the steps statistics of simulations (via Eq. and 
is in this graph indicated by the dashed red line. Although 
the errors caused by our model (Dh + 0^c) i s oru y a ^ ew 
percent, it is still important to realize that step correlation 
are negligible when the actual diffusion constants (solid dots) 
meet the dashed line. 
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C. Finite size fractal sphere packings: Levy Glass 

So far, we have investigated unbounded holey systems. 
In practice, one often deals with finite size system where 
transport can be even more complicated to understand. 
For homogenous media it is, for example, well known that 
the transport in slabs deviates from diffusion when sam- 
ple thickness L is less than about 10 times the transport 
mean free path [46-48 . For heterogeneous systems, we 
anticipate that additional complications arise when the 
sizes of non-scattering regions are on the order of the sam- 
ple size. This section therefore aims at discussing the re- 
lation between transport in unbounded media, as treated 
in previous sections, and transport in finite size media. 
We will show that the theory of diffusion constants is not 
directly applicable to bounded systems when the size of 
heterogeneities is on the order of system dimensions. 

This issue of bounded heterogeneous systems is partic- 
ularly relevant to the ongoing discussion on anomalous 
transport in Levy Glass [3ffl36l 149] , Levy Glasses are 
turbid materials with strong (fractal) spatial heterogene- 
ity in the density of scatterers. The heterogeneity is en- 
gineered and controlled by the embedding non-scattering 
regions that follow a power law size distribution into a 
turbid medium. More specifically, the non-scattering re- 
gions are constituted by glass spheres of sizes that are 
(ideally) exponentially sampled from some smallest ra- 
dius r m i n to an upper radius r max . Assuming one sphere 
crossing per step and a negligible contribution from the 
inter-sphere media, setting the number density of spheres 
n ~ r -(«+2) should produce a step length distribution 
that fall off as £ _ ( Q; + 1 ) [49 . On scales on the order of the 
maximal sphere size, as in slabs of thickness L = 2r max , 
the transport is then believed to be close to that of a 
Levy walk. Although the probabilistic theory presented 
here represents a significant contribution to the topic of 
Levy Glass design, this falls outside the main focus of 
this article. A view on this specific but important mat- 
ter is, nonetheless, given in Appendix |B] Here, we will 
instead stick to the more general question of the relation 
between transport in unbounded and bounded media. 

The system studied in the previous section is, in fact, 
an unbounded Levy Glass with a = 1. Its composition 
closely resembles the Levy Glass systems that have been 
studied experimentally [31j [35] and simulation-wise [36] . 
Our study has shown that step correlations play an im- 
portant role for transport in the unbounded systems, and 
that the resulting diffusion constant for the system was 
about 2111 |Ltm 2 /ps. Looking back at Fig. u\ and the 
evolution of the mean square displacement (MSD), it ap- 
pears that the onset of diffusion occurs surprisingly early. 
The MSD starts to grow approximately linear with time 
already after 2 ps. At this time, the MSD is about 0.027 
mm 2 , meaning that the average walker displacement is 
(roughly) on the order of 150 |im. From this, it may 
be tempting to conclude that transport is diffusive af- 
ter some 150 |im, and that the transmission through a 
Levy Glass of thickness L = 2r max = 400 |Ltm indeed 



will be diffusive. In fact, Groth et al. [36 have claimed 
that transport through Levy Glass follows regular diffu- 
sion. To check whether the above reasoning is adequate 
or not, and to investigate the claim of Groth et al., we 
have performed simulations of the transport through a 
bounded version of this sphere packing. 




FIG. 9. 2D projection of a simulated trajectory through a 
3D Levy Glass slab. In this case, the random walker crossed 
one of the largest spheres. To obtained average transmission 
properties, walkers are injected randomly on top of the slab 
(the slab consists of 9507676 spheres, the thickness being L = 
2r max + e = 402 ujn, cj) = 0.7005). Analysis of the sphere 
packing reveals that the exponential spacing model applies 
equally well as in the unbounded case, and that t s = A = 
5.89 [im indeed renders E[N] very close to one. 

Fig. [9] illustrates our simulations of transport through 
bounded systems. Our conclusion from these simulations 
is that transport through Levy Glass cannot be described 
by diffusion. Along with a significant amount of ballistic 
and quasi-ballistic transmission, we observe major devi- 
ations from diffusion theory also in mean transmission 
time and long-time decay constant (lifetime). Diffusion 
theory of light transport is well established, and in the 
diffusive regime it is, for example, well known [50] that 
the mean transmission time in follows 



(L + 2z e ) 2 
2D 



and that the decay time constant tjj is given by 



TD 



(L + 2z e f 

7T 2 D 



(25) 



(26) 



In the above formulas, L denotes slab thickness, and 
z e = the so called extrapolation length. If transport 
through the slab is diffusive with D = 2111 |Ltm 2 /ps, 
macroscopic transport should be identical to a random 
walk of independent and exponential distributed steps 
with average length i t = 3D/v « 33 |Ltm. We would thus 
expected a mean transmission time of about ijj = 15 ps 
and a decay constant of about td = 9 ps. In contrast, 
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as shown in Fig. [lOj quenched simulations result in a 
mean transmission time of 17.4 ps and decay constant of 
12.5 ps. In fact, the observed combination of t and r is 
in fact not consistent with any diffusion constant. Fur- 
thermore, the significant amount of ballistic and quasi- 
ballistic light alone indicates that diffusion cannot well 
capture the transport. We therefore conclude that trans- 
port in Levy Glass is not governed by diffusion. Instead, 
we found that the dynamics is well reproduced by, what 
we call, a quasi- annealed simulation in which we mimic 
the step distribution but effectively remove step anti- 
correlations. In this straight-forward simulation, being 
a numerical correspondent to our probabilistic theory, 
steps are generated on the basis of the exponential spac- 
ing model under the requirement that when the walker is 
about to a make a chord, this chord cannot be longer than 
the distance to the slab boundary (see Appendix C]for a 
more detailed description). The step length is thus po- 
sition and direction dependent, but step correlations are 
heavily suppressed. Although it is out of scope of this 
article to investigate this interesting result in detail, it 
appears as if the presence of boundaries in the quenched 
system strongly alter the role of step correlations. 



Probability density 
0.06 

t b 

0.04 
0.02 



Quenched MC 
t = 17.29 ± 0.10 ps 
r= 12.40 ±0.02 ps 

Quasi-annealed MC 
t = 16.97 ±0.01 ps 
r= 11.62 ±0.04 ps 




10 15 20 
Time-of-flight (ps) 



FIG. 10. Transmission dynamics for a Levy Glass, i.e. a 
bounded fractal sphere packing (slab thickness L = 402 ]±m, 
maximal sphere diameter 400 \im). Important characteristics 
such as mean transmission time t and decay constant r show 
that the dynamics is not consistent with diffusion (see main 
text). This conclusion is also supported by the significant 
amount of quasi-ballistic transmission (truly ballistic compo- 
nent at tb = L/v is not part of the shown distribution). The 
dynamics is instead well reproduced by a quasi- annealed sim- 
ulation based on the exponential spacing model. This mutual 
agreement indicates that step correlations play a less impor- 
tant role here than in the unbounded case. Statistics are 
based on 3 and 10 repetitions of 10 6 walkers for the quenched 
and quasi-annealed simulations, respectively (mean times and 
decay constant are stated in the graph as the mean ± stan- 
dard error). The transmission were in both cases about 12%, 
while the the ballistic fraction of the transmission was 0.6% 
and 0.15% (higher in the quenched case). 



IV. CONCLUSIONS 

We have presented a probabilistic theory for transport 
in systems with quenched heterogeneous distribution of 
scatterers. The focus lied on random systems, and via a 
simple model of the quenched disorder we have been able 
to derive analytical expressions for, notably, the mean 
step, the mean squared step, the diffusion constant and 
void crossing probabilities. Here we repeat two key re- 
sults. First, we have shown that the mean step is given 
by (Eq. § 



E[S] 



1 



4 



and thus equals the homogenized (exponential) scatter- 
ing mean free path being dependent only of void fill- 
ing fraction 0, not their size and shape. Second, for holey 
systems in the form of random sphere packings we have 
also shown (in the limit of weak step correlations) that 
the asymptotic diffusion constant is (Eq. [23]) 



D 



1 



-vih 



• 6 x -v x — —-^ 
V 6 E[(] 



In cases where step correlations are non-negligible, un- 
der the reasonable assumption that steps are not posi- 
tively correlated, this expression can function as an up- 
per bound on the diffusion constant. All our theoretical 
results are well supported by Monte Carlo simulations of 
random walks in holey systems based on random sphere 
packings. We have also shown that the relation between 
transport in unbounded and bounded system are par- 
ticularly complex when heterogeneities are on the order 
the system size (especially regarding the onset of regular 
diffusion and the role of step correlations). 

Given the abundance of turbid materials where scat- 
terers are not homogeneously distributed, we believe that 
both our viewpoint and our results can be of value for a 
wide range of topics - from fundamental work on light 
transport to applied spectroscopy of heterogeneous me- 
dia such as biological tissues, food products, and powder 
compacts. In particular, viewing transport as a quenched 
holey random walk may turn out to be an important 
approach for work directed towards the understanding 
of the optics and spectroscopy of complex porous me- 
dia with broad and/or anisotropic pore size distribution, 
such as pharmaceutical tablets [STJ [52] . In this context, 
it should be noted that generalization of our work to 
systems with impedance mismatch and anisotropic het- 
erogeneity should be rather straight forward. Account 
for refractive index mismatch is, for example, standard in 
the study of light transport [53, 54 . Finally, as discussed 
earlier, a theoretical account of quenched systems with 
strong step correlations is more challenging, but equally 
important. 
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Appendix A: Derivations for E[S 2 ] 



In this appendix, we provide a derivation of the expres- 
sion for the mean squared step given in the main text. 
We start by expanding the mean squared step into three 
components: 



E[S 2 



E [^turbid] 



^[^void] 



■2E[Si 



turbid ^void I 



(Al) 



Since ^[£ t 2 ur bid] = ^1 P er definition, the challenge is 
the last two terms. As mentioned in the main text, we 
will make the important assumption that the length of 
the different chords involved in a step is independent on 
the number of spheres crossed (this means, e.g., that we 
assume that the spacing between two spheres are inde- 
pendent of the sizes of these spheres). By splitting the 
outcome of SVoid into cases of different number of sphere 
crossings, we can calculate the two terms from the rules 
of total expectation values. Starting with E[S^ oid ], we 
have 



^ v 2 oid] = E P » xi? [ 5 voidl^ = n] 

n=0 

oo n 
= Y / PnXE[(Y / Q) 2 ] = 

oo 

= J2PnX(nE[C 2 }+n(n-l)E[C} 2 ) 
= E[N]E[( 2 } + (e[N 2 } - E[N]JE[C] 



<A2) 



When calculating i£[Sturbid5'void] 5 we must remember 
that they are dependent variables. We can break this de- 
pendence by, again, splitting into cases of different num- 



E [turbid SVoid] = ^] Pji x E [turbid ^void l^ci 



i=0 

oo 



Y,Pn x nEfc] x £[S turbid |iV = n] 

AL 



i=0 

oo 



Y,Pn xn£[C] x (n+1) 



i=0 



A + L 



E[C]^rx{E[N 2 }+E[N} 



A + i s 



2E[(] 



2E[(] 



Ms 
A + L 



AL (P. £ 



E[N] 2 + E[N] 



A + i s \A 2 A 



= 2E[Q^,^L + A) 



= 2E[(] x 



A + L A 2 



(A3) 



The step to reach the third line, involving a calculation of 
^[Sturbidl^V = w], deserves a comment. This term tells us 
the average inter- void path made given a certain number 
of crossings. Having n crossings is the same as knowing 
that the step was longer than a sum of n spacings, but 
that the step in the could not cross the (n + l):th spacing. 
This event can be written as 



n 

< ^turbid <(J2 A i) + A 
i=l 



(A4) 



Noting that the sum of the n spacing is an Erlang 
distribution (sum on n independent exponential distri- 
butions), we can calculate ^[^turbidl^ = n] by do- 
ing three-dimensional integration over the joint distribu- 
tion of three random variables (the Erlang distribution, 
A n+ i G Exp (A) and 5t U rbid £ Exp(^ s )). The conditional 
expectation value is reached by integrating over the space 
corr esponding to the event TV = n, as expressed in Eq. 
A4 (note that the joint density function must be renor- 
malized by the probability of the event). The result is, 
in essence, the Gamma function integral. A more elegant 
way to solve the problem is, however, to see the problem 
as an encounter of spacings A^, where the memoryless 
nature of the random process makes it possible to treat 
each and every crossing independently. If a walker man- 
ages to cross a certain spacing, the best estimate of the 
spacing crossed is i£[A|St U rbid > A]. On the other hand, 
when the walker finally does not manage to cross a spac- 
ing, the step made in this last part (after the last void 
crossing) is on average E'[S'turbid|A > ^turbid]- These two 
terms can be calculated by means of a two-dimensional 
integral over their joint probability density function (the 
symmetry of the problem makes one term follow from the 
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other). 

E[A\S turUd > A] 



Jff5xp(£,S)d£dS 
JJJ P (£,S)£d5 
AL 



A + £ s 

= ^[^turbid|A > turbid] (A5) 

So, when a walker manages to cross a inter-sphere spac- 
ing, that spacing is, on average, i£[A|Sturbid > A]. As 
emphasized above, the memoryless nature of exponen- 
tial distribution allows us to forget the total inter-sphere 
step already made, and make the same conclusion at af- 
ter each crossing. When the walker finally fails to cross a 
spacing, the average step taken after the last void cross- 
ing was [^turbid | A > Sturbid]- Given that there was n 
crossings in total, the expected path is therefore 

E[S turhid \N = n]=nx E[A\S turhid > A] 
+ E[S\A > S turhid ] 
AL 



= (n + l)- 



£, 



(A6) 



We reach the expression Eq. 20 by summing the three 
contributions and reducing an appearing quadratic poly- 
nomial 



E[S 2 } = E[S* urhid ] 
= 2£ 2 S + <t>i h 



^[^void] 
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E[(] 
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E[C] 



f 2 



4£[C] 



P. 
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Appendix B: Levy Glass design 

When Levy Glass were introduced by Barthelemy et 
al. [31 , the scatterer concentration was selected aim- 
ing at achieving, on average, one scattering event be- 
tween sphere crossings (i.e. E[N] = 1, where N states 
the number of sphere crossings in a step). The design 
principle was not elaborated in more detail, and exactly 
how filling fraction and scattering strength affect cross- 
ings statistics was not known. The inter-sphere transport 
mean free path was in fact chosen to match the smallest 
sphere diameter, i.e. £ t = 2r m i n = 5 |im. The proba- 
bilistic theory presented in this work allows a selection 
of the scattering strength that more strictly fulfills the 
ideal Levy Glass design criteria. From Eq. [T7| we see 
that E[N] = 1 is reached when £ s = A. As can be seen 



in Fig. |6j this result is confirmed by our simulations. 
In this respect (using Eq. 19 to express A), the ideal 
scattering strength is therefore 



£- 



ideal 



E[(] 



(Bl) 



Returning to the experimental work in Ref. [31] , the fill- 
ing fraction, as can be calculated from recipe given in the 
supplemental information, was about 67%. Via calcula- 
tion of the mean chord, we find that the ideal £ s for the 
samples used to investigate transmission scaling ranges 
from 9.5 |Ltm (L = 550 |im, L = 2r max being the Levy 
Glass slab thickness) to 4.2 |Ltm (L = 50 |Ltm). Although 
the utilized scattering strength clearly is on the right or- 
der of magnitude, setting £ s = 2r m i n is not universally 
ideal. That the ideal £ s varies with Levy Glass thick- 
ness is related to the fact that thicker samples include 
the use of larger spheres that fill space more efficiently 
than smaller ones. If the filling fraction <j> is kept con- 
stant, the average sphere spacing will therefore increase 
with thickness. It is thus important to realize that stud- 
ies of a series of Levy Glass with varying L = 2r^ but 
with constant <ft and £ s comes with an inherent mismatch 
with the ideal Levy Glass design principle. The impli- 
cation this finding has on the design of scaling exper- 
iments remains to be elaborated. For completeness, it 
should also be noted that scattering is not fully isotropic 
in the discussed experiments. The titania nanoparticles 
used gives a anisotropy factor of g « 0.6, meaning that 
£ s = (1 — g)£ t 7^ £ t . This is an additional aspect that 
should be considered in the future. 

In a more recent experimental paper, Burresi et al. [35] 
reported on weak localization of light (coherent backscat- 
tering) in Levy glass. There, samples had a filling frac- 
tion of about 70% and were manufactured using spheres 
of sizes from r m i n = 2.5 |Ltm up to r max = 115 |Ltm 
(L « 230 |xm). Our theory shows that the ideal £ s 
is around 6.3 |Ltm and 3.8 |Ltm for the studied a = 1 
and a = 1.5 samples, respectively. In the article, refer- 
ence experiments where scatterers were dispersed homo- 
geneously rendered £ t = 19 |Ltm, indicating that the inter- 
sphere scattering strength was l t = 19/(1 — 0) = 5.7 pin. 
The experiments done were thus in close agreement with 
the ideal design principle (again, disregarding that scat- 
tering was anisotropic with g = 0.6). 

Finally, a recent paper by Groth, Akhmerov and 
Beenakker [36] investigated transmission scaling in Levy 
Glasses via simulations of random walks in sphere pack- 
ings. There, simulations are stated to run with £ s = 
^min/2 and it is said that this choice was made to ensure 
that "there is, on average, one scattering event between 
leaving and entering a sphere" (i.e., to ensure E[N = 1]). 
Given the findings of our present work, we believe that 
this level of scattering is too strong to render E[N] = 1 
(at least for relevant filling fractions). Above, we re- 
ported that experimental Levy Glass systems have had 
a filling fraction of about 70%, and that the ideal scat- 
tering mean free path always is significantly larger than 
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the radius of the smallest sphere. Assuming that Groth 
et al. are not studying very different systems, we believe 
that using t s = r m i n /2 brings them far from the ideal 
situation where E[N] = 1. However, since utilized filling 
fractions were not stated, we cannot be fully quantitative. 
Nonetheless, it is important to note that while setting i s 
very small will take us close to the desired power 
law decay in the step length distribution (multi crossing 
getting increasing unlikely, i.e., P(N > 1) negligible), 
such a procedure will induce strong step correlations that 
affect transport. 



Appendix C: The quasi-annealed random walk 

The proposed quasi- annealed random walk aims at 
mimicking the step length distribution of the quenched 
system while removing step correlations related to the 
quenched disorder. After random walker initiation, and 
after each scattering event, a random step Sturbid £ 
Exp (-£ s ) to be travelled through the turbid component 
is generated. The length of this step is first compared 
to the distance to the boundary, along the current 
direction of propagation. If Sturbid > the walker 
leaves the sample and a transmission or reflection time 
is registered. If, on the other hand, Sturbid < 5turbid 
is instead compared to a randomly generated exponen- 
tial distributed distance to the "closest virtual sphere", 
A s G Exp(A) (cf. the exponential spacing model de- 
scribed in Sect. 



IID). If Si 



turbid 



< A s , scattering will 
take place, and the procedure will start over when a 
new direction and a new step Sturbid nave been gen- 
erated. If 5turbid > A s , the walker will cross a non- 
scattering region, after which the remaining part of the 
step 5turbid — will be used as the new step forward 
(to be compared with the updated A5 and a new sphere 
spacing). The length travelled through the void (a ran- 
dom chord (") is generated based on the distribution of 
chords predicted by equilibrium considerations, i.e. the 



density function given in Eq. [l2j The related cumulative 
distribution functions (CDF) is 



F c (x) = P(( < x) = 1 - P(C > x) 
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Solving y = F{(x), we find 



x = F-\y) 
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(C5) 



We use this inverse CDF to generate the random chord 
length C, setting 



(C6) 



where U is a random variable uniformly distributed in 
[0,1]. Note that a generated chord is only accepted if 
it is smaller than the current distance to the boundary. 
If a generated chord is larger than the distance to the 
boundary it is discarded, and the generation is remade. 
This means that long steps to some extent are under- 
represented compared to the equilibrium consideration 
that lies behind the chord length distribution. At this 
stage, it cannot be ruled out that this has some impact on 
quantities such as mean transit time and decay rate (life- 
time). The role of step correlations in transport through 
quenched disorder certainly deserves further attention in 
the future. 



[1] Thomas, G. & Stamnes, K. Radiative Transfer in the 
Atmosphere and Ocean (Cambridge University Press, 
2002). 

[2] Mishchenko, M., Travis, L. & Lacis, A. Multiple scatter- 
ing of light by particles: radiative transfer and coherent 
backs cattering (Cambridge University Press, 2006). 

[3] Davis, A. B. & Marshak, A. Solar radiation transport in 
the cloudy atmosphere: a 3D perspective on observations 
and climate impacts. Rep. Prog. Phys. 73, 026801 (2010). 

[4] Tuchin, V. V. Tissue Optics: Light Scatter- 
ing Methods and Instruments for Medical Diagnosis 
(SPIE/International Society for Optical Engineering, 
2007). 

[5] Welch, A. & Gemert, M. Optical- Thermal Response of 
Laser- Irradiated Tissue (Springer, 2010), 2 edn. 



[6] Wang, L. & Wu, H. Biomedical Optics: Principles and 
Imaging (Wiley, 2007). 

[7] Berne, B. & Pecora, R. Dynamic Light Scattering: With 
Applications to Chemistry, Biology, and Physics (Dover 
Publications, 2000). 

[8] Reich, G. Near- infrared spectroscopy and imaging: Basic 
principles and pharmaceutical applications. Adv. Drug 
Deliver. Rev. 57, 1109-1143 (2005). 

[9] Siesler, H., Ozaki, Y., Kawata, S. & Heise, H. Near- 
Infrared Spectroscopy: Principles, Instruments, Applica- 
tions (John Wiley & Sons, 2008). 
[10] Shi, Z. Q. & Anderson, C. A. Pharmaceutical applica- 
tions of separation of absorption and scattering in near- 
infrared spectroscopy (NIRS). J. Pharm. Sci. 99, 4766- 
4783 (2010). 



13 



[11] Shaw, R. A., Kostinski, A. B. & Lanterman, D. D. Super- 
exponential extinction of radiation in a negatively cor- 
related random medium. J. Quant. Spectrosc. Radiat. 
Transfer 75, 13-20 (2002). 

[12] Davis, A. B. & Marshak, A. Photon propagation in het- 
erogeneous optical media with spatial correlations: en- 
hanced mean- free-paths and wider-than-exponential free- 
path distributions. J. Quant. Spectrosc. Radiat. Transfer 
84, 3-34 (2004). 

[13] Davis, A. B. & Mineev-Weinstein, M. B. Radiation prop- 
agation in random media: From positive to negative cor- 
relations in high-frequency fluctuations. J. Quant. Spec- 
trosc. Radiat. Transfer 112, 632-645 (2011). 

[14] Bal, G. & Jing, W. Fluctuation theory for radiative 
transfer in random media. J. Quant. Spectrosc. Radiat. 
Transfer 112, 660-670 (2011). 

[15] Lovejoy, S., Davis, A., Gabriel, P., Schertzer, D. & 
Austin, G. L. Discrete angle radiative transfer 1. scal- 
ing and similarity, universality and diffusion. J. Geophys. 
Res. 95, 11699-11715 (1990). 

[16] Davis, A., Marshak, A. & Pfeilsticker, K. Anomalous 
Levy photon diffusion theory: toward a new parameter- 
ization of shortwave transport in cloudy columns. 9th 
ARM Science Team Meeting Proceedings (1999). 

[17] Scholl, T., Pfeilsticker, K., Davis, A. B., Klein Baltink, 
H., Crewell, S., Lhnert, U., Simmer, C., Meywerk, J. 
& Quante, M. Path length distributions for solar pho- 
tons under cloudy skies: Comparison of measured first 
and second moments with predictions from classical and 
anomalous diffusion theories. J. Geophys. Res. Ill, 
D12211 (2006). 

[18] Liu, H., Chance, B., Hielscher, A., Jacques, S. & Tit- 
tel, F. Influence of blood-vessels on the measurement of 
hemoglobin oxygenation as determined by time-resolved 
reflectance spectroscopy. Med. Phys. 22, 1209-1217 
(1995). 

[19] Firbank, M., Arridge, S. R., Schweiger, M. & Delpy, 
D. T. An investigation of light transport through scatter- 
ing bodies with non-scattering regions. Phys. Med. Biol. 
41, 767-783 (1996). 

[20] Hielscher, A., Alcouffe, R. & Barbour, R. Comparison of 
finite-difference transport and diffusion calculations for 
photon migration in homogeneous and heterogeneous tis- 
sues. Phys. Med. Biol. 43, 1285-1302 (1998). 

[21] Boas, D., Culver, J., Stott, J. & Dunn, A. Three dimen- 
sional Monte Carlo code for photon migration through 
complex heterogeneous media including the adult human 
head. Opt. Express 10, 159-170 (2002). 

[22] Bal, G. Particle transport through scattering regions 
with clear layers and inclusions. J. Comput. Phys. 180, 
659-685 (2002). 

[23] Shah, N., Cerussi, A. E., Jakubowski, D., Hsiang, D., 
Butler, J. & Tromberg, B. Spatial variations in optical 
and, physiological properties of healthy breast tissue. J. 
Biomed. Opt. 9, 534-540 (2004). 

[24] Gibson, A., Hebden, J. & Arridge, S. Recent advances 
in diffuse optical imaging. Phys. Med. Biol. 50, R1-R43 
(2005). 

[25] Ntziachristos, V., Ripoll, J., Wang, L. H. V. & 
Weissleder, R. Looking and listening to light: the evo- 
lution of whole-body photonic imaging. Nat. Biotechnol. 
23, 313-320 (2005). 

[26] Svensson, T., Andersson-Engels, S., Einarsdottfr, M. & 
Svanberg, K. In vivo optical characterization of human 



prostate tissue using near-infrared time-resolved spec- 
troscopy. J. Biomed. Opt. 12, 014022 (2007). 

[27] Dogariu, A., Uozumi, J. & Asakura, T. Enhancement 
of the backscattered intensity from fractal aggregates. 
Waves Random Media 2, 259-263 (1992). 

[28] Sorensen, CM. Light scattering by fractal aggregates: 
A review. Aerosol Sci. Technol. 35, 648-687 (2001). 

[29] Durian, D. J., Weitz, D. A. & Pine, D. J. Multiple light- 
scattering probes of foam structure and dynamics. Sci- 
ence 252, 686-688 (1991). 

[30] Gittings, A. S., Bandyopadhyay, R. & Durian, D. J. Pho- 
ton channelling in foams. Europhys. Lett. 65, 414-419 
(2004). 

[31] Barthelemy, P., Bertolotti, J. & Wiersma, D. S. A Levy 
flight for light. Nature 453, 495-498 (2008). 

[32] Burioni, R., Caniparoli, L. & Vezzani, A. Levy walks and 
scaling in quenched disordered media. Phys. Rev. E 81, 
060101 (2010). 

[33] Barthelemy, P., Bertolotti, J., Vynck, K., Lepri, S. & 
Wiersma, D. S. Role of quenching on superdiffusive trans- 
port in two-dimensional random media. Phys. Rev. E 82, 
011101 (2010). 

[34] Buonsante, P., Burioni, R. & Vezzani, A. Transport 
and scaling in quenched two- and three-dimensional Levy 
quasicrystals. Phys. Rev. E 84, 021105 (2011). 

[35] Burresi, M., Radhalakshmi, V., Savo, R., Bertolotti, J., 
Vynck, K. & Wiersma, D. S. Weak localization of light 
in superdiffusive random systems. Phys. Rev. Lett. 108, 
110604 (2012). 

[36] Groth, C. W., Akhmerov, A. R. & Beenakker, C. W. J. 
Transmission probability through a Levy glass and com- 
parison with a levy walks. Phys. Rev. E 85, 021138 
(2012). 

[37] Lagendijk, A. & van Tiggelen, B. A. Resonant multiple 
scattering of light. Phys. Rep. 270, 143-215 (1996). 

[38] van Rossum, M. C. W. & Nieuwenhuizen, T. M. Multiple 
scattering of classical waves: microscopy, mesoscopy, and 
diffusion. Rev. Mod. Phys. 71, 313-371 (1999). 

[39] Bouchaud, J. -P. & Georges, A. Anomalous diffusion in 
disordered media: Statistical mechanisms, models and 
physical applications. Phys. Rep. 195, 127-293 (1990). 

[40] Ben-Avraham, D. & Havlin, S. Diffusion and reactions 
in fractals and disordered systems (Cambridge University 
Press, 2000). 

[41] Olson, G. L., Miller, D. S., Larsen, E. W. & Morel, J. E. 

Chord length distributions in binary stochastic media in 

two and three dimensions. J. Quant. Spectrosc. Radiat. 

Transfer 101, 269-283 (2006). 
[42] De Kruijf, W. & Kloosterman, J. On the average chord 

length in reactor physics. Ann. Nucl. Energy 30, 549-553 

(2003). 

[43] Torquato, S. & Lu, B. Chord- length distribution function 
for two-phase random media. Phys. Rev. E 47, 2950- 
2953 (1993). 

[44] Gueron, D. & Mazzolo, A. Properties of chord length 
distributions across ordered and disordered packing of 
hard disks. Phys. Rev. E 68, 066117 (2003). 

[45] Torquato, S. Random Heterogeneous Materials: Mi- 
crostructure and Macroscopic Properties (Springer, 
2001). 

[46] Yoo, K., Liu, F. & Alfano, R. When does the 
diffusion-approximation fail to describe photon trans- 
port in random-media? Phys. Rev. Lett. 64, 2647-2650 
(1990). 



14 



[47] Kop, R. H. J., de Vries, P., Sprik, R. & Lagendijk, A. 
Observation of anomalous transport of strongly multiple 
scattered light in thin disordered slabs. Phys. Rev. Lett. 
79, 4369-4372 (1997). 

[48] Elaloufi, R., Carminati, R. & Greffet, J.- J. Diffusive-to- 
ballistic transition in dynamic light transmission through 
thin scattering slabs: a radiative transfer approach. J. 
Opt. Soc. Am. A 21, 1430-1437 (2004). 

[49] Bertolotti, J., Vynck, K., Pattelli, L., Barthelemy, P., 
Lepri, S. & Wiersma, D. S. Engineering disorder in su- 
perdiffusive Levy glasses. Adv. Funct. Mater. 20, 965-968 
(2010). 

[50] Vellekoop, I. M., Lodahl, P. & Lagendijk, A. Determina- 
tion of the diffusion constant using phase-sensitive mea- 



surements. Phys. Rev. E 71, 056604 (2005). 

[51] Svensson, T., Alerstam, E., Johansson, J. & Andersson- 
Engels, S. Optical porosimetry and investigations of the 
porosity experienced by light interacting with porous me- 
dia. Opt. Lett. 35, 1740-1742 (2010). 

[52] Alerstam, E. & Svensson, T. Observation of anisotropic 
diffusion of light in compacted granular porous materials. 
Phys. Rev. E 85, 040301 (2012). 

[53] Wang, L., Jacques, S. & Zheng, L. MCML Monte Carlo 
modeling of light transport in multilayered tissues. Corn- 
put. Meth. Prog. Bio. 47, 131-146 (1995). 

[54] Sadjadi, Z. & Miri, M. Diffusive transport of light in two- 
dimensional granular materials. Phys. Rev. E 84, 051305 
(2011). 



